Regresión Lineal Múltiple

TP Integrador





Alumnos

Nicolás Dominutti, Carlos Suarez Gurruchaga, Hernán Telechea









set.seed(101)

COVID = read.csv("COVID.txt", sep="")
dim(COVID)
## [1] 139  35
summary(COVID)
##     geoId             CntrName            muertes           casos          
##  Length:139         Length:139         Min.   :     0   Min.   :      2.0  
##  Class :character   Class :character   1st Qu.:    12   1st Qu.:    725.5  
##  Mode  :character   Mode  :character   Median :    66   Median :   3333.0  
##                                        Mean   :  2518   Mean   :  37621.0  
##                                        3rd Qu.:   548   3rd Qu.:  20625.0  
##                                        Max.   :100442   Max.   :1699933.0  
##  muertes.permil     casos.permil       l10muertes.permil   Hombres80      
##  Min.   :  0.000   Min.   :    0.949   Min.   :0.0000    Min.   :0.08737  
##  1st Qu.:  2.214   1st Qu.:   70.054   1st Qu.:0.5062    1st Qu.:0.37722  
##  Median :  6.074   Median :  277.583   Median :0.8497    Median :1.07801  
##  Mean   : 52.735   Mean   : 1133.633   Mean   :1.0282    Mean   :1.59580  
##  3rd Qu.: 33.787   3rd Qu.: 1596.139   3rd Qu.:1.5414    3rd Qu.:2.57442  
##  Max.   :819.817   Max.   :17596.220   Max.   :2.9142    Max.   :6.16625  
##    Mujeres80          Pobla80          Pobla65          PoblaMid    
##  Min.   : 0.1639   Min.   :0.1308   Min.   : 1.085   Min.   :47.42  
##  1st Qu.: 0.5983   1st Qu.:0.4885   1st Qu.: 3.284   1st Qu.:60.15  
##  Median : 1.4421   Median :1.2710   Median : 6.450   Median :64.70  
##  Mean   : 2.6441   Mean   :2.1200   Mean   : 9.223   Mean   :63.58  
##  3rd Qu.: 4.5155   3rd Qu.:3.7069   3rd Qu.:14.920   3rd Qu.:66.90  
##  Max.   :10.5564   Max.   :8.3613   Max.   :27.576   Max.   :85.09  
##    PoblaData           PoblaDens           Mujeres          Urbano      
##  Min.   : 0.001102   Min.   : 0.02041   Min.   :24.50   Min.   : 13.03  
##  1st Qu.: 0.047427   1st Qu.: 0.33375   1st Qu.:49.72   1st Qu.: 44.62  
##  Median : 0.111232   Median : 0.80729   Median :50.40   Median : 62.45  
##  Mean   : 0.485381   Mean   : 1.87767   Mean   :49.89   Mean   : 60.83  
##  3rd Qu.: 0.333277   3rd Qu.: 1.46399   3rd Qu.:51.01   3rd Qu.: 78.83  
##  Max.   :13.927300   Max.   :79.52998   Max.   :54.54   Max.   :100.00  
##    ExpectVida      NeontlMort       DisMort          Lesion      
##  Min.   :49.84   Min.   : 0.90   Min.   : 7.80   Min.   : 2.600  
##  1st Qu.:64.57   1st Qu.: 3.40   1st Qu.:14.30   1st Qu.: 5.650  
##  Median :71.11   Median : 8.70   Median :17.90   Median : 8.700  
##  Mean   :70.06   Mean   :12.30   Mean   :18.30   Mean   : 8.827  
##  3rd Qu.:75.14   3rd Qu.:20.45   3rd Qu.:22.35   3rd Qu.:10.700  
##  Max.   :81.70   Max.   :42.00   Max.   :30.60   Max.   :28.400  
##    EnfNoTrans       EnfTrans          PBI            Tuberculosis  
##  Min.   :26.00   Min.   : 1.30   Min.   :  0.7568   Min.   :  1.0  
##  1st Qu.:48.50   1st Qu.: 5.20   1st Qu.:  4.1778   1st Qu.: 13.0  
##  Median :74.80   Median :11.20   Median : 12.3515   Median : 45.0  
##  Mean   :69.39   Mean   :21.78   Mean   : 20.0941   Mean   :106.5  
##  3rd Qu.:88.75   3rd Qu.:38.00   3rd Qu.: 28.3214   3rd Qu.:150.0  
##  Max.   :95.20   Max.   :65.30   Max.   :123.2139   Max.   :611.0  
##     Diabetes         Medicos            Camas         ImmunSaramp   
##  Min.   : 1.000   Min.   :0.01767   Min.   : 0.200   Min.   :30.00  
##  1st Qu.: 5.050   1st Qu.:0.28709   1st Qu.: 1.147   1st Qu.:85.00  
##  Median : 6.500   Median :1.38303   Median : 2.171   Median :93.00  
##  Mean   : 7.425   Mean   :1.65391   Mean   : 2.980   Mean   :88.21  
##  3rd Qu.: 9.200   3rd Qu.:2.77551   3rd Qu.: 3.984   3rd Qu.:97.00  
##  Max.   :22.000   Max.   :5.71011   Max.   :13.796   Max.   :99.00  
##    TempMarzo         HipTen.H        HipTen.M         HipTen     
##  Min.   :-18.72   Min.   : 9.00   Min.   :12.40   Min.   :10.70  
##  1st Qu.:  4.93   1st Qu.:25.85   1st Qu.:25.85   1st Qu.:26.07  
##  Median : 19.99   Median :32.60   Median :32.00   Median :32.00  
##  Mean   : 15.38   Mean   :32.80   Mean   :31.98   Mean   :32.39  
##  3rd Qu.: 25.25   3rd Qu.:40.40   3rd Qu.:38.00   3rd Qu.:39.83  
##  Max.   : 30.63   Max.   :54.60   Max.   :50.60   Max.   :50.60  
##       BCG             BCGf               Tiempo      
##  Min.   :0.0000   Length:139         Min.   :  0.00  
##  1st Qu.:1.0000   Class :character   1st Qu.: 57.00  
##  Median :1.0000   Mode  :character   Median : 65.00  
##  Mean   :0.9496                      Mean   : 63.15  
##  3rd Qu.:1.0000                      3rd Qu.: 75.00  
##  Max.   :1.0000                      Max.   :136.00

Hay 35 variables analizadas para las 139 filas de datos. Una de las variables es la lista de países. 3 de las variables son cualitativas, mientras que las demás son cuantitativas.

attach(COVID)

par(mfrow=c(1, 2))

boxplot(muertes.permil ~ BCGf, 
        main = "Muertes por millón de habitantes",
        ylab = "Muertes.permil", col=c("thistle","wheat"), 
        cex = 0.7, xlab = "Estrategia de inmunización",
        names = c("selectiva", "todos"))

boxplot(l10muertes.permil ~ BCGf, 
        main = "Muertes por millón [log]",
        ylab = "log10muertes.permil", col=c("thistle","wheat"),
        cex = 0.7, xlab = "Estrategia de inmunización", 
        names = c("selectiva", "todos"))

Ante esta comparación de boxplots, conviene considerar la variables l10muertes.permil porque la relación entre BCGf y muertes.permil no parece ser muy lineal y estaría considerando muchos valores como potencialmente atípicos cuando en realidad no lo son. Al aplicar logaritmo, apenas un dato es considerado extraño, mientras que el resto de los datos parecen estar más claramente distribuidos. Me doy cuenta mucho más fácilmente que tener una estrategia de vacunación contra BCG selectiva podría generar un mayor número de muertes en las poblaciones.


Continuamos el análisis con la variable l10muertes.permil.

#cor(COVID[,-c(1,2,34)])

Las relaciones más fuertes entre variables son: Pobla80 con Hombre80 y Mujer80, que tiene mucho sentido dado que la población es una función de los hombres y mujeres mayores a 80 años. El total de mujeres en la población total es mayor a la de los hombres, por eso su correlación es ~0.99 (la de los hombres es ~0.98). En cambio, las muertes no parecen tener relación alguna con otras enfermedades o la densidad poblacional (correlaciones menores al módulo de 0.05).

Planteamos algunas preguntas para identificar relaciones entre variables (muchas de ellas dentro de las más consultadas durante la pandemia):

plot(PBI,l10muertes.permil)

cor(PBI,l10muertes.permil)
## [1] 0.5385521

Es poco clara la relación. Si bien un alto porcentaje de los países con bajo PBI tienen un bajo total de muertes logarítmicas por millón, es poco claro el incremento (o decrecimiento) de esta variable a mayores niveles de PBI. Es débil esta relación: la correlación es positiva en ~0.53.

plot(PoblaMid,l10muertes.permil)

plot(Pobla65,l10muertes.permil)

plot(Pobla80,l10muertes.permil)

Al comparar 3 tipos de población, identificamos que, cuanto mayor es la edad promedio de la población, más relación lineal hay con respecto al total de muertes. Esto se comprueba con las correlaciones:

cor(PoblaMid,l10muertes.permil)
## [1] 0.3623718
cor(Pobla65,l10muertes.permil)
## [1] 0.6579062
cor(Pobla80,l10muertes.permil)
## [1] 0.6676371

Las 3 correlaciones son positivas pero la mayor se da con Pobla80 (0.6676).

plot(TempMarzo,l10muertes.permil)

cor(TempMarzo,l10muertes.permil)
## [1] -0.5107478

Al parecer, temporaturas más elevadas parecen generar menores chances de mortalidad, quizás porque los contagios son menores y/o el virus resiste menos. De todas formas, no es fuerte la relación (correlación es negativa en ~0.51).

mean(l10muertes.permil[TempMarzo > 20] < 1)
## [1] 0.7971014

El ~80% de los países cuya temperatura promedio en marzo fue mayor a 20 grados tuvieron menos de 1 muerte logarítmica por mil.

plot(Tiempo,TempMarzo)

mean(Tiempo[TempMarzo > 20] > 60)
## [1] 0.8405797

Sin embargo, estos datos solo son de marzo de 2020, época en la que la enfermedad aún no estaba tan propagada en países del hemisferio sur (más cálido en esa época). Por ende, en lugar de decir que hubo menos muertes porque el ambiente era más cálido, podríamos inferir que esto fue así porque la enfermedad tardó más en llegar. Con esto, también visualizamos que los países con mayores temperaturas en marzo son también los que más tardaron en registrar contagios, en promedio. De todas maneras, no sabemos si esos contagios tardaron en llegar por la temperatura o por alguna otra causa que desconocemos.

plot(Camas,l10muertes.permil)

cor(Camas,l10muertes.permil)
## [1] 0.3636445
plot(Medicos,l10muertes.permil)

cor(Medicos,l10muertes.permil)
## [1] 0.6191814
df = COVID[ ,c(5:7, 10, 11, 12, 23, 24, 25, 26, 27, 35)]
pairs(df[ ,c("l10muertes.permil", "Pobla80", "Pobla65", 
             "PBI", "Tuberculosis", "Medicos", "Camas")])

Esta relación parece contradictoria. Si bien no es fuerte, mucho no dice que tener camas ayuda a que la gente no se muera. No es un factor tan determinante. La relación entre médicos y muertes tiene un poco más de sentido. En los mismos lugares en los que mayor número de muertes hay, mayor es la presencia de médicos, quizás por el tamaño poblacional. La correlación es positiva en ~0.62.

plot(Urbano,l10muertes.permil)

cor(Urbano,l10muertes.permil)
## [1] 0.5717382

Es bastante claro que ambientes más urbanizados tienen mayor número de muertes. Esto es porque la población total también es mayor.

COVID$Incidencia = round(COVID$muertes / COVID$casos,2)
hist(COVID$Incidencia)

Vemos que, en la mayoría de los países, la incidencia fue de entre 0 y 0.05. Veamos cuáles son los que quedan hacia la cola derecha:

COVID %>% 
  filter((COVID$Incidencia>.1)) %>%
  select('CntrName')

Como punto en común, podemos decir que todos estos países pertenecen al hemisferio norte y, salvo México, Belize y Yemen, +70% de ellos son europeos.

df = COVID %>%
     mutate(alto = case_when(Incidencia > .1 ~ 1, Incidencia <= .1 ~ 0))

pairs(df %>% 
         filter(df$alto==1) %>%
         select('Pobla80','muertes.permil',
                'casos.permil','l10muertes.permil'),
         col = "darkblue")

En este pairplot, podemos ver que, dentro de los países donde la mortalidad del virus sobre los casos totales fue mayor al .1%, existe una posible relación marcada entre la proporción de habitantes con +80 años y las muertes por mil.

aux = df %>% filter(df$alto==1)

plot(y=aux$Pobla80, x=aux$muertes.permil, main='Efecto de poblaciones longevas\n sobre muertes en países con letalidad alta')
abline(h=4.5, col='red')

En este deepdive, podemos ver que se podría segmentar estos países en 2 clusters, teniendo en cuenta el marcado salto que se produce en las muertes por mil cuando la variable Pobla80 supera el 4.4%.

Tomamos otras variables que, intuitivamente, parecieran tener relación directa en muertes o, en algunos casos, en los casos de covid.

¿Existirá alguna relación detectable?

pairs(df %>%
        select('ExpectVida','Pobla80','l10muertes.permil'), col = "darkblue")

Vemos que, en países donde la expectativa de vida es mayor, se produjeron más muertes. Esto se podría explicar por la relación directa (logarítmica) que tiene esta variable con la proporción de personas mayores a 80 años presentes en este tipo de países y, cómo vimos en el punto anterior, esta variable se encuentra relacionada con las muertes.

Ahora, intentemos detectar estas relaciones lineales de las que hablamos previamente, calculando el índice de correlación de Pearson:

corrplot(
  cor(select_if(df, is.numeric)),
  type='lower',
  #addCoef.col="black",
  number.cex=0.75,
  mar=c(0,0,1,0),
  title = 'Matriz de correlaciones'
  )


lm.covid = lm(l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
                                  Tuberculosis + Camas + TempMarzo + PBI, 
              data = COVID)

lm.covid_2 = lm(muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
                                  Tuberculosis + Camas + TempMarzo + PBI, 
              data = COVID)

summary(lm.covid)
## 
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
##     Tuberculosis + Camas + TempMarzo + PBI, data = COVID)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.08088 -0.29012 -0.03314  0.27929  1.29267 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.5981046  0.2012630   2.972  0.00352 ** 
## PoblaDens    -0.0111873  0.0063386  -1.765  0.07990 .  
## Pobla80       0.2040913  0.0338369   6.032 1.55e-08 ***
## Urbano        0.0082702  0.0025767   3.210  0.00167 ** 
## Tuberculosis -0.0005815  0.0003431  -1.695  0.09246 .  
## Camas        -0.1101291  0.0263087  -4.186 5.17e-05 ***
## TempMarzo    -0.0143117  0.0057068  -2.508  0.01337 *  
## PBI           0.0062466  0.0027181   2.298  0.02314 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4765 on 131 degrees of freedom
## Multiple R-squared:  0.6239, Adjusted R-squared:  0.6038 
## F-statistic: 31.04 on 7 and 131 DF,  p-value: < 2.2e-16
summary(lm.covid_2)
## 
## Call:
## lm(formula = muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
##     Tuberculosis + Camas + TempMarzo + PBI, data = COVID)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -230.52  -42.58   -8.56   21.49  650.26 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   10.84325   43.17327   0.251 0.802087    
## PoblaDens     -1.03060    1.35971  -0.758 0.449840    
## Pobla80       40.79237    7.25840   5.620  1.1e-07 ***
## Urbano         0.30373    0.55274   0.550 0.583594    
## Tuberculosis   0.05273    0.07359   0.717 0.474932    
## Camas        -22.21224    5.64354  -3.936 0.000134 ***
## TempMarzo     -1.59048    1.22418  -1.299 0.196149    
## PBI            1.18989    0.58307   2.041 0.043285 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 102.2 on 131 degrees of freedom
## Multiple R-squared:  0.3765, Adjusted R-squared:  0.3432 
## F-statistic:  11.3 on 7 and 131 DF,  p-value: 3.772e-11

Optamos por utilizar como variable target l10muertes.permil, debido a su mayor R2 y cantidad de regresoras significativas para el modelo.

leverage_sorted=sort(hatvalues(lm.covid), decreasing = T)

boxplot(leverage_sorted, 
        main = "LEVERAGE DE LAS OBSERVACIONES PARA DISTINTOS PAISES",
        col = "blue")

leverage_sorted[1:10]
##        SG        QA        JP        MN        CA        BY        LU        BN 
## 0.9451793 0.2946915 0.2162310 0.1745735 0.1486524 0.1356584 0.1347664 0.1336935 
##        GA        LS 
## 0.1333837 0.1307780

Se observa que el modelo le asigna un leverage muy elevado a la observacion que representa a Singapur, seguida por las que representan a Qatar, Japan, Mongolia, Canada, Bielorusia, Luxemburgo, Brunei_Darussalam, Gabon y Lesotho, respectivamente.

plot(lm.covid, 4)

En esta representancion de las distancias de Cook, observamos tambien que las observaciones de Singapur, Japan y Qatar, presentan elevadas distancias de Cook, lo cual reafirma que tienen un leverage elevado, ya que la misma es proporcional a el leverage de la observacion, asi como tambien al residuo de la misma. Por lo que ahora queda pendiente analizar los residuos de estas observaciones para descartar o no que sean posibles outliers, o simplemente observaciones muy influyentes para nuestro modelo.

plot(lm.covid, 3)

Al analizar el aspecto de los residuos del modelo, pareciera haber un grado de HETEROCEDASTICIDAD de los residuos, pero no visualizamos ninguna estructura en los mismos.

cor.test(abs(lm.covid$residuals), lm.covid$fitted.values, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  abs(lm.covid$residuals) and lm.covid$fitted.values
## S = 367626, p-value = 0.0355
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1786362

Con un pvalue pequeño rechazamos H0, y por tanto las varianzas de los residuos son distintas, hay HETEROCEDASTICIDAD.

plot(lm.covid, 2)

ks.test(lm.covid$residuals, "pnorm",
                             mean(lm.covid$residuals), 
                             sd(lm.covid$residuals))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  lm.covid$residuals
## D = 0.048945, p-value = 0.8931
## alternative hypothesis: two-sided

Con un pvalue de 0.89, nuestros residuos tienen una distribucion que se asemeja a una normal.

plot(lm.covid, 5)

Al analizar los residuos en funcion del Leverage de las observaciones, vemos que la observacion correspondiente a Singapur, si bien tiene un leverage elevado como mencionamos en el item anterior, no posee un residuo alto, y por tanto (por el momento) preferimos dejarla en el modelo. Tambien vemos resaltadas las observaciones Qatar y Japan, las cuales tambien poseen un alto valor de Leverage, pero su residuo no parece ser alto y por tanto no serian outliers.

step(lm.covid, direction = "both")
## Start:  AIC=-198.35
## l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + Tuberculosis + 
##     Camas + TempMarzo + PBI
## 
##                Df Sum of Sq    RSS     AIC
## <none>                      29.738 -198.34
## - Tuberculosis  1    0.6522 30.390 -197.33
## - PoblaDens     1    0.7071 30.445 -197.08
## - PBI           1    1.1989 30.937 -194.85
## - TempMarzo     1    1.4277 31.166 -193.83
## - Urbano        1    2.3385 32.076 -189.82
## - Camas         1    3.9778 33.716 -182.90
## - Pobla80       1    8.2586 37.996 -166.28
## 
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
##     Tuberculosis + Camas + TempMarzo + PBI, data = COVID)
## 
## Coefficients:
##  (Intercept)     PoblaDens       Pobla80        Urbano  Tuberculosis  
##    0.5981046    -0.0111873     0.2040913     0.0082702    -0.0005815  
##        Camas     TempMarzo           PBI  
##   -0.1101291    -0.0143117     0.0062466

Realizamos un analisis de los residuos del modelo en funcion de las covariables utilizadas a partir del criterio “STEPWISE” y llegamos a la conclusion que el mejor modelo es el que utiliza todas las variables propuestas en el item 3), sin quitar ninguna.

Respecto a las covariables que mas influencia cobran para predecir el target, se encuentran el promedio de poblacion masculina y femenina mayor a 80 años en primer medida, seguido por la cantidad de camas de hospital cada 1000 personas.

ANALISIS DE OUTLIERS POR REGRESION ROBUSTA

Vamos a plantear una regresion robusta para chequear si los betas estimados por OLS dan parecidos, o hay observaciones que nos estan afectando la regresion y no las encontramos en los analisis hechos.

control = lmrobdet.control(bb = 0.5, efficiency = 0.85, family = "bisquare")



modelo_MM_covid = lmrobdetMM(l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
                                  Tuberculosis + Camas + TempMarzo + PBI, 
                             data = COVID, control = control)

summary(modelo_MM_covid)
## 
## Call:
## lmrobdetMM(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + Tuberculosis + 
##     Camas + TempMarzo + PBI, data = COVID, control = control)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.24220 -0.27260 -0.03026  0.31517  1.31445 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.3849847  0.2222787   1.732  0.08563 .  
## PoblaDens    -0.0125395  0.0063513  -1.974  0.05045 .  
## Pobla80       0.2524079  0.0373247   6.762 4.07e-10 ***
## Urbano        0.0091258  0.0027644   3.301  0.00124 ** 
## Tuberculosis -0.0004969  0.0003576  -1.390  0.16697    
## Camas        -0.1163727  0.0286581  -4.061 8.37e-05 ***
## TempMarzo    -0.0091945  0.0063197  -1.455  0.14809    
## PBI           0.0060534  0.0029396   2.059  0.04145 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 0.4764 
## Multiple R-squared:  0.6361, Adjusted R-squared:  0.6166 
## Convergence in 25 IRWLS iterations
(lm.covid$coefficients)
##   (Intercept)     PoblaDens       Pobla80        Urbano  Tuberculosis 
##  0.5981046383 -0.0111873307  0.2040912643  0.0082702004 -0.0005814871 
##         Camas     TempMarzo           PBI 
## -0.1101291038 -0.0143117235  0.0062466456
(modelo_MM_covid$coefficients)
##   (Intercept)     PoblaDens       Pobla80        Urbano  Tuberculosis 
##  0.3849846701 -0.0125395234  0.2524078877  0.0091257877 -0.0004969057 
##         Camas     TempMarzo           PBI 
## -0.1163726595 -0.0091945121  0.0060533946

Los beta estimados por ambos metodos son muy semejantes, lo que nos da la pauta de que NO tenemos outliers que esten viciando nuestro modelo original de LS, respecto a los coeficientes que no dan estadisiticamente significativos para el metodo robusto, puede deberse a que el mismo por el metodo que emplea para su calculo, sus intervalos son mas largos que los que devuelve OLS y por tanto, puede ocurrir que el 0 caiga dentro del intervalo, dandonos a creer que una variable no tiene informacion util para darnos, cuando no es asi. Por esto vamos a continuar con el metodo de OLS, dado que no tenemos outliers y es el metodo que nos devuelve los intervalos de los betas mas pequeños.

plot(y = modelo_MM_covid$rweights, 
     x = modelo_MM_covid$residuals, col = "blue",
     ylim = c(0, 1), xlab = "RESIDUOS", ylab = "RWEIGHTS",
     main = "RWEIGHTS vs RESIDUOS MODELO ROBUSTO")

abline(h = 0, col = "red", lwd = 2)

Reafirmamos que no tenemos observaciones que sean outliers, ya que el metodo MMestimadores, no tiene ninguna observacion que le asigne un peso de 0

df2 = COVID[, c("PoblaDens","Pobla80", "Urbano", "Tuberculosis", "Camas", "TempMarzo", "PBI")]


pairs(df2)

Probamos crear una interaccion entre las variables PBI y Camas

df2[, "PBIcamas"] = df2$Camas * df2$PBI

lm.covid_interaccion = lm(l10muertes.permil ~ PoblaDens + Pobla80 + 
                            Urbano + Tuberculosis + Camas + 
                            TempMarzo + PBI + PBIcamas, data = df2)

summary(lm.covid_interaccion)
## 
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
##     Tuberculosis + Camas + TempMarzo + PBI + PBIcamas, data = df2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.07242 -0.29205 -0.03691  0.28372  1.29236 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.6201428  0.2148048   2.887  0.00456 ** 
## PoblaDens    -0.0111062  0.0063665  -1.744  0.08344 .  
## Pobla80       0.2004483  0.0360438   5.561 1.46e-07 ***
## Urbano        0.0084492  0.0026531   3.185  0.00181 ** 
## Tuberculosis -0.0005942  0.0003468  -1.713  0.08907 .  
## Camas        -0.1181358  0.0374609  -3.154  0.00200 ** 
## TempMarzo    -0.0147272  0.0058904  -2.500  0.01366 *  
## PBI           0.0051970  0.0044247   1.175  0.24233    
## PBIcamas      0.0003444  0.0011433   0.301  0.76369    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4781 on 130 degrees of freedom
## Multiple R-squared:  0.6241, Adjusted R-squared:  0.601 
## F-statistic: 26.98 on 8 and 130 DF,  p-value: < 2.2e-16

NO observamos una mejoria del R2 respecto al modelo original.

df2[,"PBIpobla80"] = df2$PBI * df2$Pobla80

lm.covid_interaccion_2 = lm(l10muertes.permil ~ PoblaDens + Pobla80 + 
                            Urbano + Tuberculosis + Camas + 
                            TempMarzo + PBI + PBIpobla80, data = df2)

summary(lm.covid_interaccion_2)
## 
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
##     Tuberculosis + Camas + TempMarzo + PBI + PBIpobla80, data = df2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.08353 -0.28711 -0.05352  0.31835  1.30154 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.6286418  0.2019852   3.112 0.002283 ** 
## PoblaDens    -0.0118262  0.0063384  -1.866 0.064322 .  
## Pobla80       0.1436949  0.0565796   2.540 0.012272 *  
## Urbano        0.0089970  0.0026267   3.425 0.000823 ***
## Tuberculosis -0.0006721  0.0003488  -1.927 0.056176 .  
## Camas        -0.1048197  0.0265341  -3.950 0.000127 ***
## TempMarzo    -0.0139004  0.0056986  -2.439 0.016066 *  
## PBI           0.0031349  0.0035807   0.875 0.382926    
## PBIpobla80    0.0016874  0.0012690   1.330 0.185937    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4751 on 130 degrees of freedom
## Multiple R-squared:  0.6289, Adjusted R-squared:  0.6061 
## F-statistic: 27.54 on 8 and 130 DF,  p-value: < 2.2e-16

Tampoco observamos un aumento significativo en el R2 al interaccionar PBI con Pobla80.

Continuamos probando interacciones en busca de una mejora del R2.

df2[, "PBItemp"] = df2$TempMarzo * df2$PBI

lm.covid_interaccion_3 = lm(l10muertes.permil ~ PoblaDens + Pobla80 + 
                            Urbano + Tuberculosis + Camas + 
                            TempMarzo + PBI + PBItemp, data = df2)

summary(lm.covid_interaccion_3)
## 
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
##     Tuberculosis + Camas + TempMarzo + PBI + PBItemp, data = df2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.03547 -0.28698 -0.06073  0.30877  1.30492 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.4306473  0.2233762   1.928 0.056048 .  
## PoblaDens    -0.0075220  0.0066628  -1.129 0.260995    
## Pobla80       0.1820208  0.0360832   5.044  1.5e-06 ***
## Urbano        0.0089768  0.0025934   3.461 0.000728 ***
## Tuberculosis -0.0006461  0.0003429  -1.884 0.061764 .  
## Camas        -0.0994583  0.0268896  -3.699 0.000319 ***
## TempMarzo    -0.0067999  0.0072204  -0.942 0.348057    
## PBI           0.0113280  0.0040551   2.793 0.006003 ** 
## PBItemp      -0.0003461  0.0002061  -1.679 0.095517 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4732 on 130 degrees of freedom
## Multiple R-squared:  0.6319, Adjusted R-squared:  0.6092 
## F-statistic: 27.89 on 8 and 130 DF,  p-value: < 2.2e-16
df2[, "PBIurbano"] = df2$Urbano * df2$PBI

lm.covid_interaccion_4 = lm(l10muertes.permil ~ PoblaDens + Pobla80 + 
                            Urbano + Tuberculosis + Camas + 
                            TempMarzo + PBI + PBIurbano, data = df2)

summary(lm.covid_interaccion_4)
## 
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
##     Tuberculosis + Camas + TempMarzo + PBI + PBIurbano, data = df2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.05951 -0.28139 -0.02248  0.31277  1.27442 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.4882083  0.2176972   2.243 0.026616 *  
## PoblaDens    -0.0086414  0.0066163  -1.306 0.193831    
## Pobla80       0.1931678  0.0347706   5.555 1.50e-07 ***
## Urbano        0.0096886  0.0027905   3.472 0.000702 ***
## Tuberculosis -0.0005501  0.0003430  -1.604 0.111162    
## Camas        -0.1138959  0.0263971  -4.315 3.14e-05 ***
## TempMarzo    -0.0140060  0.0056964  -2.459 0.015258 *  
## PBI           0.0230117  0.0131379   1.752 0.082208 .  
## PBIurbano    -0.0001902  0.0001458  -1.304 0.194486    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4752 on 130 degrees of freedom
## Multiple R-squared:  0.6287, Adjusted R-squared:  0.6059 
## F-statistic: 27.52 on 8 and 130 DF,  p-value: < 2.2e-16
df2[, "PBI^2"] = df2$PBI * df2$PBI

lm.covid_interaccion_5 = lm(l10muertes.permil ~ PoblaDens + Pobla80 + 
                            Urbano + Tuberculosis + Camas + 
                            TempMarzo + PBI + PBI^2, data = df2)

summary(lm.covid_interaccion_5)
## 
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
##     Tuberculosis + Camas + TempMarzo + PBI + PBI^2, data = df2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.08088 -0.29012 -0.03314  0.27929  1.29267 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.5981046  0.2012630   2.972  0.00352 ** 
## PoblaDens    -0.0111873  0.0063386  -1.765  0.07990 .  
## Pobla80       0.2040913  0.0338369   6.032 1.55e-08 ***
## Urbano        0.0082702  0.0025767   3.210  0.00167 ** 
## Tuberculosis -0.0005815  0.0003431  -1.695  0.09246 .  
## Camas        -0.1101291  0.0263087  -4.186 5.17e-05 ***
## TempMarzo    -0.0143117  0.0057068  -2.508  0.01337 *  
## PBI           0.0062466  0.0027181   2.298  0.02314 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4765 on 131 degrees of freedom
## Multiple R-squared:  0.6239, Adjusted R-squared:  0.6038 
## F-statistic: 31.04 on 7 and 131 DF,  p-value: < 2.2e-16

Si bien no conseguimos una mejora muy significativa, la mejor interaccion que conseguimos fue entre PBI y Urbano. Vamos a aplicar nuevamente el criterio de seleccion STEPWISE, para ver si adoptamos esta interaccion o mantenemos nuestro modelo origil.

step(lm.covid_interaccion_4, direction = "both")
## Start:  AIC=-198.15
## l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + Tuberculosis + 
##     Camas + TempMarzo + PBI + PBIurbano
## 
##                Df Sum of Sq    RSS     AIC
## - PBIurbano     1    0.3840 29.738 -198.34
## - PoblaDens     1    0.3852 29.739 -198.34
## <none>                      29.354 -198.15
## - Tuberculosis  1    0.5809 29.935 -197.43
## - PBI           1    0.6927 30.047 -196.91
## - TempMarzo     1    1.3650 30.719 -193.83
## - Urbano        1    2.7219 32.076 -187.83
## - Camas         1    4.2036 33.557 -181.55
## - Pobla80       1    6.9689 36.323 -170.54
## 
## Step:  AIC=-198.35
## l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + Tuberculosis + 
##     Camas + TempMarzo + PBI
## 
##                Df Sum of Sq    RSS     AIC
## <none>                      29.738 -198.34
## + PBIurbano     1    0.3840 29.354 -198.15
## - Tuberculosis  1    0.6522 30.390 -197.33
## - PoblaDens     1    0.7071 30.445 -197.08
## - PBI           1    1.1989 30.937 -194.85
## - TempMarzo     1    1.4277 31.166 -193.83
## - Urbano        1    2.3385 32.076 -189.82
## - Camas         1    3.9778 33.716 -182.90
## - Pobla80       1    8.2586 37.996 -166.28
## 
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
##     Tuberculosis + Camas + TempMarzo + PBI, data = df2)
## 
## Coefficients:
##  (Intercept)     PoblaDens       Pobla80        Urbano  Tuberculosis  
##    0.5981046    -0.0111873     0.2040913     0.0082702    -0.0005815  
##        Camas     TempMarzo           PBI  
##   -0.1101291    -0.0143117     0.0062466

Nuestro analisis de residuos por STEPWISE considera que es mejor no tomar en cuenta la interaccion, asi que vamos a continuar con nuestro modelo original.

int_prediccion = predict(lm.covid, newdata = COVID, interval = c("predict"))

plot(y = lm.covid$fitted.values, 
     x = COVID$l10muertes.permil, 
     col = "darkblue", ylim = c(-2, 5), cex = .8,
     xlab = "VALORES OBSERVADOS", ylab = "VALORES PREDICHOS")

# Intervalos de prediccion
points(x = COVID$l10muertes.permil, 
      y = int_prediccion[, 2], col = "darkred", lwd =3, cex = 0.5)

points(x = COVID$l10muertes.permil, 
       y = int_prediccion[, 3], col = "darkred", lwd =3, cex = 0.5)

legend("topleft", c("INTERVALO DE PREDICCION", "OBSERVACION"), fill = c("darkred", "darkblue"))


PROBLEMAS DE COLINEALIDAD:

Por construccion sabemos que geoID nos otorga la misma informacion que CntrName.

Por construccion sabemos que casos y casos.permil, estan altamente correlacionadas.

Por construccion sabemos que muertes y muertes.permil, estan altamente correlacionadas.

Por construccion sabemos que Pobla80, esta relacionado altamente con Mujeres80 y Hombres80, ya que es el promedio de los dos anteriores.

Por construccion sabemos tambien que EnfNoTrans y EnfTrans guardan una colinealidad por construccion.

Por construccion sabemos tambien que HipTen esta relacionada altamente con HipTen.H e HipTen.M, ya que es el promedio de los dos anteriores.

Por construccion sabemos tambien que BCG y BCGf tienen correlacion 1.

Sobre estos casos, una posible solucion seria optar por alguna de las variables que estan correlacionadas entre si, eliminando las demas, que no aportan informacion util nueva. Por lo tanto vamos a adoptar quedarnos con CntrName, casos.permil, muertes.permil, Pobla80, EnfTrans, HipTen, BCG, eliminando las demas.

A su vez, tambien vamos a realizar un analisis de multicolinealidad a partir del Factor de inflacion de la varianza (VIF)

df_for_vif = COVID[, -c(1, 2, 3, 4, 5, 8, 9, 21, 30, 31, 34, 36)]
lm_for_vif = lm(l10muertes.permil ~., data = df_for_vif)

barplot(vif(lm_for_vif)[vif(lm_for_vif)>5], main = "VIF DE LAS COVARIABLES",
                                             col = "blue", cex.names = 0.65)
abline(h = 5, lwd = 2, col = "red")

Recordemos que el VIF es una medida de inflación de los coeficientes de la variable independiente debido a multicolinealidad entre las variables regresoras. Un VIF que excede 5 requiere investigación de multicolinealidad y un VIF de >10 indica multicolinealidad.

Investigaremos luego los siguientes casos de alto VIF más elevados:

cov_altoVIF = names(vif(lm_for_vif)[vif(lm_for_vif)>5])

data.frame(t(cor(df_for_vif$Pobla80, df_for_vif[, cov_altoVIF], 
               method = "spearman")))
data.frame(t(cor(df_for_vif$Pobla65, df_for_vif[, cov_altoVIF], 
               method = "spearman")))
data.frame(t(cor(df_for_vif$ExpectVida, df_for_vif[, cov_altoVIF], 
               method = "spearman")))
data.frame(t(cor(df_for_vif$EnfTrans, df_for_vif[, cov_altoVIF], 
               method = "spearman")))
df3 = COVID[, -c(1, 3, 4, 8, 9, 21, 30, 31, 34, 36)]

La covariable CntrName, podriamos codificarla como una variable dummie para poder alimentar el modelo, pero optamos por no incluirla en la regresion, ya que realizar ese proceso, aumentaria aun mas la dimensionalidad de nuestro problema que ya de por si es alta, y relacion costo beneficio, no creemos que lo valga.

covariables_l10muertes=df3[, -c(1, 2)]
covariables_muertes_permil=df3[, -c(1, 4)]

modelo_muertes_permil = lm(muertes.permil~. , data = covariables_muertes_permil)
summary(modelo_muertes_permil)
## 
## Call:
## lm(formula = muertes.permil ~ ., data = covariables_muertes_permil)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -225.55  -26.13   -5.17   20.92  354.33 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -14.32051  607.74146  -0.024  0.98124    
## casos.permil    0.02909    0.00578   5.033 1.80e-06 ***
## Pobla80        63.45564   22.45560   2.826  0.00556 ** 
## Pobla65       -15.27977    7.39576  -2.066  0.04107 *  
## PoblaMid       -2.86408    3.15996  -0.906  0.36664    
## PoblaData       1.12118    4.86447   0.230  0.81813    
## PoblaDens      -2.08770    1.27482  -1.638  0.10423    
## Mujeres         7.12569    4.30409   1.656  0.10054    
## Urbano         -0.05985    0.53843  -0.111  0.91168    
## ExpectVida      2.53633    4.25542   0.596  0.55233    
## NeontlMort      1.06751    1.88626   0.566  0.57254    
## DisMort         0.05534    3.12626   0.018  0.98591    
## Lesion         -3.35844    2.71100  -1.239  0.21793    
## EnfTrans       -1.21433    1.88601  -0.644  0.52095    
## PBI             0.07827    0.74049   0.106  0.91600    
## Tuberculosis    0.02925    0.08207   0.356  0.72217    
## Diabetes       -2.95779    2.61345  -1.132  0.26009    
## Medicos       -10.35879   12.50693  -0.828  0.40925    
## Camas          -5.09684    5.41314  -0.942  0.34839    
## ImmunSaramp     0.18782    0.74213   0.253  0.80066    
## TempMarzo       0.85855    1.27358   0.674  0.50158    
## HipTen         -0.14586    1.23836  -0.118  0.90645    
## BCG          -209.21052   37.58929  -5.566 1.73e-07 ***
## Tiempo         -0.60811    0.50497  -1.204  0.23096    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 81.8 on 115 degrees of freedom
## Multiple R-squared:  0.6494, Adjusted R-squared:  0.5793 
## F-statistic: 9.261 on 23 and 115 DF,  p-value: < 2.2e-16
modelo_l10muertes = lm(l10muertes.permil~. , data = covariables_l10muertes)
summary(modelo_l10muertes)
## 
## Call:
## lm(formula = l10muertes.permil ~ ., data = covariables_l10muertes)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.52890 -0.22590  0.00776  0.25171  0.98278 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -4.015e+00  3.230e+00  -1.243  0.21634    
## casos.permil  1.670e-04  3.071e-05   5.436  3.1e-07 ***
## Pobla80       6.175e-02  1.193e-01   0.517  0.60584    
## Pobla65       2.951e-02  3.930e-02   0.751  0.45423    
## PoblaMid      1.632e-02  1.679e-02   0.972  0.33323    
## PoblaData     4.407e-03  2.585e-02   0.170  0.86492    
## PoblaDens    -1.879e-02  6.774e-03  -2.773  0.00647 ** 
## Mujeres       5.051e-02  2.287e-02   2.208  0.02922 *  
## Urbano        6.648e-03  2.861e-03   2.323  0.02192 *  
## ExpectVida    2.247e-02  2.261e-02   0.994  0.32244    
## NeontlMort    5.602e-03  1.002e-02   0.559  0.57736    
## DisMort       5.238e-04  1.661e-02   0.032  0.97490    
## Lesion       -2.224e-03  1.441e-02  -0.154  0.87761    
## EnfTrans      7.238e-03  1.002e-02   0.722  0.47164    
## PBI          -2.326e-03  3.935e-03  -0.591  0.55567    
## Tuberculosis -7.972e-04  4.361e-04  -1.828  0.07016 .  
## Diabetes      1.209e-02  1.389e-02   0.871  0.38575    
## Medicos       2.294e-02  6.646e-02   0.345  0.73062    
## Camas        -6.068e-02  2.877e-02  -2.109  0.03708 *  
## ImmunSaramp  -7.843e-03  3.944e-03  -1.989  0.04911 *  
## TempMarzo    -7.626e-03  6.768e-03  -1.127  0.26217    
## HipTen       -1.524e-03  6.581e-03  -0.232  0.81731    
## BCG          -2.400e-01  1.998e-01  -1.201  0.23204    
## Tiempo        4.089e-05  2.683e-03   0.015  0.98787    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4347 on 115 degrees of freedom
## Multiple R-squared:  0.7252, Adjusted R-squared:  0.6702 
## F-statistic: 13.19 on 23 and 115 DF,  p-value: < 2.2e-16

Observamos que utilizando la transformacion logaritmica del target obtenemos el R2 mas elevado, asi como tambien mas estimadores estadisticamente distintos de cero. Aun asi, observamos que debido a la gran varianza que tienen nuestros estimadores, hay muchos coeficientes que no logran ser significativos, esto puede deberse a los problemas de colinealidad que mencionamos previamente y aun no eliminamos.

Mejor modelo hasta el momento: modelo_l10muertes.


Vamos a aplicar el metodo de seleccion de variables stepwise, utilizando como metrica a minimizar el criterio de AKAIKE, para intentar disminuir los problemas de colinealidad y con ello mejorar nuestro modelo.

wise_modelo_l10muertes = step(modelo_l10muertes, direction = c("both"), 
                                                 trace = F)

summary(wise_modelo_l10muertes)
## 
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla80 + PoblaDens + 
##     Mujeres + Urbano + ExpectVida + Tuberculosis + Camas + ImmunSaramp + 
##     BCG, data = covariables_l10muertes)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.57911 -0.26035 -0.00065  0.26120  0.91573 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.292e+00  1.095e+00  -2.093  0.03831 *  
## casos.permil  1.604e-04  2.546e-05   6.299 4.42e-09 ***
## Pobla80       1.373e-01  3.729e-02   3.682  0.00034 ***
## PoblaDens    -1.811e-02  5.472e-03  -3.310  0.00121 ** 
## Mujeres       4.485e-02  1.522e-02   2.946  0.00383 ** 
## Urbano        5.304e-03  2.370e-03   2.238  0.02695 *  
## ExpectVida    2.121e-02  1.056e-02   2.009  0.04666 *  
## Tuberculosis -5.820e-04  3.573e-04  -1.629  0.10585    
## Camas        -4.556e-02  2.185e-02  -2.086  0.03900 *  
## ImmunSaramp  -7.726e-03  3.489e-03  -2.214  0.02858 *  
## BCG          -3.009e-01  1.791e-01  -1.680  0.09544 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4219 on 128 degrees of freedom
## Multiple R-squared:  0.7118, Adjusted R-squared:  0.6893 
## F-statistic: 31.61 on 10 and 128 DF,  p-value: < 2.2e-16
barplot(vif(wise_modelo_l10muertes), main = "VIF DE LAS COVARIABLES",
                                     col = "blue", cex.names = 0.65, 
                                     ylim=c(0,6), las=2)
abline(h = 5, lwd = 2, col = "red")

Disminuimos la cantidad de covariables de nuestro modelo de 23 a 10, el R2 aumento un poco, el RSE disminuyo tambien y, al parecer, eliminamos el problema de multicolinealidad. A su vez vemos tambien que el numero de betas estimados que son significativamente distintos de cero aumento considerablemente, con lo que tenemos un mejor modelo, con menor varianza en sus coeficientes.

wise_modelo_muertes_permil = step(modelo_muertes_permil, direction = c("both"),
                                                         trace = F)

summary(wise_modelo_muertes_permil)
## 
## Call:
## lm(formula = muertes.permil ~ casos.permil + Pobla80 + Pobla65 + 
##     PoblaDens + Mujeres + Camas + BCG + Tiempo, data = covariables_muertes_permil)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -231.74  -20.62   -1.09   19.98  372.17 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.192e+02  1.388e+02  -1.579 0.116726    
## casos.permil  2.801e-02  4.673e-03   5.994 1.89e-08 ***
## Pobla80       6.426e+01  1.776e+01   3.618 0.000423 ***
## Pobla65      -1.414e+01  5.495e+00  -2.574 0.011176 *  
## PoblaDens    -1.809e+00  1.032e+00  -1.753 0.081909 .  
## Mujeres       1.026e+01  2.761e+00   3.717 0.000298 ***
## Camas        -9.202e+00  4.020e+00  -2.289 0.023698 *  
## BCG          -2.091e+02  3.357e+01  -6.227 6.07e-09 ***
## Tiempo       -7.640e-01  3.846e-01  -1.986 0.049081 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 78.64 on 130 degrees of freedom
## Multiple R-squared:  0.6337, Adjusted R-squared:  0.6111 
## F-statistic: 28.11 on 8 and 130 DF,  p-value: < 2.2e-16

Observamos nuevamente que conviene trabajar con el target convertido, pero destacamos la mejora en el R2 y RSE respecto al modelo antes de aplicar el criterio de seleccion.

Vamos a probar utilizar un enfoque backward a ver si llegamos a un

resultado distinto

back_modelo_l10muertes = step(modelo_l10muertes, direction = c("back"),
                                                 trace = F)
summary(back_modelo_l10muertes)
## 
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla80 + PoblaDens + 
##     Mujeres + Urbano + ExpectVida + Tuberculosis + Camas + ImmunSaramp + 
##     BCG, data = covariables_l10muertes)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.57911 -0.26035 -0.00065  0.26120  0.91573 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.292e+00  1.095e+00  -2.093  0.03831 *  
## casos.permil  1.604e-04  2.546e-05   6.299 4.42e-09 ***
## Pobla80       1.373e-01  3.729e-02   3.682  0.00034 ***
## PoblaDens    -1.811e-02  5.472e-03  -3.310  0.00121 ** 
## Mujeres       4.485e-02  1.522e-02   2.946  0.00383 ** 
## Urbano        5.304e-03  2.370e-03   2.238  0.02695 *  
## ExpectVida    2.121e-02  1.056e-02   2.009  0.04666 *  
## Tuberculosis -5.820e-04  3.573e-04  -1.629  0.10585    
## Camas        -4.556e-02  2.185e-02  -2.086  0.03900 *  
## ImmunSaramp  -7.726e-03  3.489e-03  -2.214  0.02858 *  
## BCG          -3.009e-01  1.791e-01  -1.680  0.09544 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4219 on 128 degrees of freedom
## Multiple R-squared:  0.7118, Adjusted R-squared:  0.6893 
## F-statistic: 31.61 on 10 and 128 DF,  p-value: < 2.2e-16

Vemos que el R2 y el RSE coinciden con el metodo stepwise, asi como los betas estimados, con lo que concluimos que el mejor modelo que podemos conseguir a este punto es wise_modelo_l10muertes.


METODOS DE REGULARIZACION

Vamos a probar ahora reducir la dimensionalidad de nuestro problema a partir de metodos de regularizacion.

xx = model.matrix(l10muertes.permil ~., data = covariables_l10muertes)[, -1]

lassocv_modelo_l10muertes <- cv.glmnet(xx, y = l10muertes.permil, 
                                         alpha = 1,
                                         nfolds = 5,
                                         family = "gaussian",
                                         type.measure = "mse")

Obtenemos los mejores valores de lambda para aplicar regularizacion lasso, usando como metrica el error cuadratico medio.

lassocv_modelo_l10muertes$lambda.min
## [1] 0.003990394
lassocv_modelo_l10muertes$lambda.1se
## [1] 0.3470641

Ahora aplicamos la regularizacion lasso con el lambda obtenido para ver que covariables sobreviven.

lasso_modelo_l10muertes_min <- glmnet(xx, y = l10muertes.permil, 
                                         alpha = 1,
                                         lambda = lassocv_modelo_l10muertes$lambda.min)

lasso_modelo_l10muertes_1sd <- glmnet(xx, y = l10muertes.permil, 
                                         alpha = 1,
                                         lambda = lassocv_modelo_l10muertes$lambda.1se)

cbind(coef(lasso_modelo_l10muertes_min),
      coef(lasso_modelo_l10muertes_1sd))
## 24 x 2 sparse Matrix of class "dgCMatrix"
##                         s0           s0
## (Intercept)  -1.7245118950 2.506694e-01
## casos.permil  0.0001601537 4.219719e-06
## Pobla80       0.0701753098 5.355739e-02
## Pobla65       0.0232764713 .           
## PoblaMid      0.0040921465 .           
## PoblaData     0.0043040149 .           
## PoblaDens    -0.0165366840 .           
## Mujeres       0.0386786518 .           
## Urbano        0.0059785151 .           
## ExpectVida    0.0125988089 9.408924e-03
## NeontlMort    0.0040718750 .           
## DisMort      -0.0039817202 .           
## Lesion       -0.0053461327 .           
## EnfTrans      .            .           
## PBI          -0.0009665838 .           
## Tuberculosis -0.0006294161 .           
## Diabetes      0.0083184355 .           
## Medicos       .            .           
## Camas        -0.0533819815 .           
## ImmunSaramp  -0.0067932573 .           
## TempMarzo    -0.0057097097 .           
## HipTen       -0.0005048359 .           
## BCG          -0.2552166366 .           
## Tiempo        .            .

Notamos que usando la regularizacion con lambda.min, no logramos reducir demasiado la dimensionalidad del problema, por lo que vamos a adoptar directamente el criterio de seleccion por lambda.1se


CRITERIO DE LAMBDA MIN + 1se

Sobreviven las covariables casos.permil, Pobla80, ExpectVida.

Si bien el metodo nos quita la variable BCG, por el tipo de problema que se esta analizando, consideramos que es importante, con lo que vamos a conservala de todos modos.

ls_covid_lasso_final_1se = lm(l10muertes.permil ~ casos.permil + Pobla80 + ExpectVida + 
                       BCG,  data = covariables_l10muertes)

summary(ls_covid_lasso_final_1se)
## 
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla80 + ExpectVida + 
##     BCG, data = covariables_l10muertes)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.78427 -0.31310  0.04202  0.29393  1.17570 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.102e-01  5.543e-01  -0.379   0.7051    
## casos.permil  1.222e-04  2.194e-05   5.571 1.34e-07 ***
## Pobla80       1.618e-01  2.751e-02   5.881 3.08e-08 ***
## ExpectVida    1.750e-02  7.897e-03   2.216   0.0284 *  
## BCG          -4.939e-01  1.912e-01  -2.584   0.0108 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4634 on 134 degrees of freedom
## Multiple R-squared:  0.6361, Adjusted R-squared:  0.6253 
## F-statistic: 58.56 on 4 and 134 DF,  p-value: < 2.2e-16

Notamos que al utilizar el metodo de regularizacion lasso, nuestro R2 NO supera al obtenido a partir del criterio stepwise. Por lo tanto vamos a conservar el modelo obtenido con stepwise, como el definitivo.

Entonces por el momento nuestro mejor modelo es el wise_modelo_l10muertes, vamos a analizar sus residuos para entender la performance

plot(wise_modelo_l10muertes)

Vemos la clara presencia de un outlier (alto leverage y alto residuo) que está dañando la regresión. El país que causa el problema es Qatar y la razón por la que es determinado como outlier es que tiene valores extremos de casos.permil, sin embargo, es uno de los países que menos muertes por mil tuvo. Tal diferencia nos llevaría a pensar que pudo haberse dado un error tanto en la recolección de los datos de casos o de muertes. También cabe destacar que es el país con mayor PBI (aunque este número pareciera ser más racional).

par(mfrow=c(1, 3))
aux = COVID[COVID$casos.permil>3000,]
barplot(aux$casos.permil, names.arg=aux$geoId, las=2, main='Casos por mil por país')
barplot(aux$muertes.permil, names.arg=aux$geoId, las=2, main='Muertes por mil por país')
barplot(aux$PBI, names.arg=aux$geoId, las=2, main='PBI por país')

Tomando este punto de evaluación como disparador, vamos a utilizar distintas técnicas para intentar mejorar la performance de nuestro mejor modelo hasta el momento.


¿Casos atípicos? Proceso de optimización del modelo


Prueba 1: drop de Qatar

Probemos ajustando una regresión sin este valor atípico, la misma debiera performar mucho mejor.

df4 = COVID[, -c(2, 3, 4, 5, 8, 9, 21, 30, 31, 34, 36)]
df4 = df4[df4$geoId!='QA',]
modelo_l10muertes_V2 = lm(l10muertes.permil~. , data = df4[,-1])
summary(modelo_l10muertes_V2)
## 
## Call:
## lm(formula = l10muertes.permil ~ ., data = df4[, -1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73444 -0.21277 -0.01593  0.20441  0.82161 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.770e-01  2.670e+00  -0.291 0.771617    
## casos.permil  3.189e-04  3.199e-05   9.969  < 2e-16 ***
## Pobla80      -3.884e-02  9.831e-02  -0.395 0.693488    
## Pobla65       6.583e-02  3.244e-02   2.030 0.044728 *  
## PoblaMid     -1.730e-03  1.391e-02  -0.124 0.901258    
## PoblaData     7.795e-03  2.111e-02   0.369 0.712627    
## PoblaDens    -2.078e-02  5.537e-03  -3.752 0.000277 ***
## Mujeres       2.887e-05  1.980e-02   0.001 0.998839    
## Urbano        3.032e-03  2.383e-03   1.272 0.205969    
## ExpectVida    1.906e-02  1.847e-02   1.032 0.304132    
## NeontlMort    6.926e-03  8.185e-03   0.846 0.399253    
## DisMort      -1.412e-04  1.356e-02  -0.010 0.991714    
## Lesion        1.494e-02  1.197e-02   1.248 0.214567    
## EnfTrans      1.734e-04  8.235e-03   0.021 0.983238    
## PBI          -5.292e-03  3.236e-03  -1.635 0.104735    
## Tuberculosis -2.979e-04  3.620e-04  -0.823 0.412187    
## Diabetes      1.592e-02  1.135e-02   1.403 0.163449    
## Medicos       4.153e-02  5.432e-02   0.765 0.446066    
## Camas        -4.001e-02  2.364e-02  -1.692 0.093315 .  
## ImmunSaramp  -6.751e-03  3.223e-03  -2.095 0.038417 *  
## TempMarzo    -1.475e-03  5.584e-03  -0.264 0.792136    
## HipTen        3.105e-03  5.407e-03   0.574 0.566871    
## BCG          -1.453e-01  1.636e-01  -0.889 0.376128    
## Tiempo        1.346e-03  2.197e-03   0.613 0.541383    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3549 on 114 degrees of freedom
## Multiple R-squared:  0.8184, Adjusted R-squared:  0.7818 
## F-statistic: 22.34 on 23 and 114 DF,  p-value: < 2.2e-16

Efectivamente logramos que la regresión aumentara su R2ajustado en +10%, pero vemos que solo 4 de los beta estimados rechazan el test-t, lo que sería un indicio probable de la presencia de multicolinealidad, vamos a correr un proceso de selección de variables sobre este modelo para intentar solucionar dicho problema.

wise_modelo_l10muertes_V2 = step(modelo_l10muertes_V2, direction = c("both"),
                                                         trace = F)

summary(wise_modelo_l10muertes_V2)
## 
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla65 + PoblaDens + 
##     Urbano + PBI + Tuberculosis + Diabetes + Camas + ImmunSaramp, 
##     data = df4[, -1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.81124 -0.19411  0.01204  0.21291  0.80697 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.6287094  0.2507573   2.507   0.0134 *  
## casos.permil  0.0003240  0.0000287  11.288  < 2e-16 ***
## Pobla65       0.0589951  0.0074604   7.908 1.06e-12 ***
## PoblaDens    -0.0232308  0.0045525  -5.103 1.18e-06 ***
## Urbano        0.0042204  0.0018903   2.233   0.0273 *  
## PBI          -0.0043099  0.0025914  -1.663   0.0987 .  
## Tuberculosis -0.0005758  0.0002664  -2.162   0.0325 *  
## Diabetes      0.0156212  0.0087189   1.792   0.0756 .  
## Camas        -0.0450035  0.0179441  -2.508   0.0134 *  
## ImmunSaramp  -0.0059518  0.0026293  -2.264   0.0253 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3431 on 128 degrees of freedom
## Multiple R-squared:  0.8094, Adjusted R-squared:  0.796 
## F-statistic:  60.4 on 9 and 128 DF,  p-value: < 2.2e-16

Vemos como al seleccionar las variables correctas pareciera que el modelo ya no tiene multicolinealidades y mejora su performance, pero nos surge una nueva duda ¿existen nuevos outliers que no hayamos detectado hasta el momento?


Prueba 2: nuevos outliers? Drop de Singapur

Como disparador en la respuesta de esta pregunta vamos inspeccionar los residuos del último modelo propuesto

plot(wise_modelo_l10muertes_V2)

Encontramos un nuevo punto con alto leverage y, pese a tener un relativamente bajo residuo estandarizado absoluta, tiene una alta distancia de Cook, por lo que investigaremos más en detalle este punto, que es Singapur (ya visto como punto extremo en el punto 3c del presente trabajo)

aux = COVID[COVID$PoblaDens>3,]
barplot(aux$PoblaDens, names.arg=aux$geoId, las=2, main='Densidad poblacional por país')

Vemos que Singapur es país con mayor densidad poblacional del set de datos, llegando estar a +79 (siendo la segunda mayor de apenas más de 12). Puede ser que, debido a este punto, nuestro modelo no pueda lograr ajustar bien a dicho país (teniendo en cuenta que el beta estimado para la variable PoblaDens es de -0.023). Vamos a probar sacando este valor, volviendo a ajustar y corriendo un proceso de selección de variables para lidiar con la multicolinealidad.

df5 = df4[df4$geoId!='SG',]
modelo_l10muertes_V3 = lm(l10muertes.permil~. , data = df5[,-1])
summary(modelo_l10muertes_V3)
## 
## Call:
## lm(formula = l10muertes.permil ~ ., data = df5[, -1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73187 -0.21917 -0.01577  0.20766  0.82485 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.408e-01  2.682e+00  -0.276   0.7829    
## casos.permil  3.191e-04  3.211e-05   9.939   <2e-16 ***
## Pobla80      -4.196e-02  9.895e-02  -0.424   0.6723    
## Pobla65       6.693e-02  3.266e-02   2.049   0.0428 *  
## PoblaMid     -1.535e-03  1.397e-02  -0.110   0.9127    
## PoblaData     6.377e-03  2.146e-02   0.297   0.7668    
## PoblaDens    -1.206e-02  2.159e-02  -0.559   0.5775    
## Mujeres       3.297e-04  1.989e-02   0.017   0.9868    
## Urbano        3.270e-03  2.459e-03   1.330   0.1863    
## ExpectVida    1.792e-02  1.874e-02   0.956   0.3409    
## NeontlMort    6.493e-03  8.280e-03   0.784   0.4346    
## DisMort      -3.072e-05  1.362e-02  -0.002   0.9982    
## Lesion        1.520e-02  1.203e-02   1.263   0.2093    
## EnfTrans      3.382e-04  8.274e-03   0.041   0.9675    
## PBI          -5.204e-03  3.255e-03  -1.599   0.1126    
## Tuberculosis -3.092e-04  3.643e-04  -0.849   0.3978    
## Diabetes      1.527e-02  1.150e-02   1.328   0.1868    
## Medicos       4.259e-02  5.457e-02   0.780   0.4368    
## Camas        -4.273e-02  2.460e-02  -1.737   0.0852 .  
## ImmunSaramp  -6.835e-03  3.241e-03  -2.109   0.0372 *  
## TempMarzo    -2.011e-03  5.749e-03  -0.350   0.7271    
## HipTen        2.999e-03  5.432e-03   0.552   0.5820    
## BCG          -1.232e-01  1.724e-01  -0.715   0.4762    
## Tiempo        1.257e-03  2.216e-03   0.567   0.5715    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3562 on 113 degrees of freedom
## Multiple R-squared:  0.8184, Adjusted R-squared:  0.7815 
## F-statistic: 22.15 on 23 and 113 DF,  p-value: < 2.2e-16
wise_modelo_muertes_permil_V3 = step(modelo_l10muertes_V3, direction = c("both"),
                                                         trace = F)
summary(wise_modelo_muertes_permil_V3)
## 
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla65 + Urbano + 
##     PBI + Tuberculosis + Diabetes + Camas + ImmunSaramp, data = df5[, 
##     -1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.81355 -0.20138  0.01161  0.21331  0.81262 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.239e-01  2.504e-01   2.492   0.0140 *  
## casos.permil  3.229e-04  2.864e-05  11.276  < 2e-16 ***
## Pobla65       5.849e-02  7.456e-03   7.844 1.49e-12 ***
## Urbano        4.662e-03  1.884e-03   2.475   0.0146 *  
## PBI          -4.229e-03  2.590e-03  -1.633   0.1050    
## Tuberculosis -5.764e-04  2.660e-04  -2.167   0.0321 *  
## Diabetes      1.328e-02  8.735e-03   1.521   0.1308    
## Camas        -4.617e-02  1.793e-02  -2.575   0.0112 *  
## ImmunSaramp  -6.255e-03  2.624e-03  -2.384   0.0186 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3426 on 128 degrees of freedom
## Multiple R-squared:  0.8097, Adjusted R-squared:  0.7978 
## F-statistic: 68.09 on 8 and 128 DF,  p-value: < 2.2e-16

Vemos que nuestro r2-ajustado mejora un poco y la gran mayoría de las variables resultan significativas

plot(wise_modelo_muertes_permil_V3)

Al ver los residuos podemos decir que ya no se ven puntos atípicos a simple vista, al parecer hemos logrado un modelo consistente y que performa bien (recordemos el salto de un modelo de 0.61 de r2-ajustado a uno de 0.7978). Pero no vamos a detenernos con lo obtenido, tenemos aun 2 approachs por probar, siendo el primero de ellos una regresión robusta


Prueba 3: regresión robusta

Para esta prueba vamos a retrotraernos al modelo antes de dropear cualquier instancia. El plan es simple, ajustar una regresión robusta y luego analizar los residuos de esta, para detectar los outliers.

control = lmrobdet.control(bb = 0.5, efficiency = 0.85, family = "bisquare")
df6 = COVID[, -c(2, 3, 4, 5, 8, 9, 21, 30, 31, 34, 36)]
modelo_MM_covid_V2 = lmrobdetMM(l10muertes.permil ~ ., 
                             data = df6[,-1], control = control)
summary(modelo_MM_covid_V2)
## 
## Call:
## lmrobdetMM(formula = l10muertes.permil ~ ., data = df6[, -1], control = control)
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.9410 -0.1828 -0.0245  0.1554  0.8794 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.748e+00  2.800e+00   0.624 0.533674    
## casos.permil  3.253e-04  3.597e-05   9.041 4.44e-15 ***
## Pobla80      -2.054e-02  1.034e-01  -0.199 0.842840    
## Pobla65       6.328e-02  3.334e-02   1.898 0.060190 .  
## PoblaMid     -9.907e-03  1.453e-02  -0.682 0.496568    
## PoblaData     3.064e-01  9.013e-02   3.400 0.000927 ***
## PoblaDens    -1.981e-02  5.600e-03  -3.537 0.000584 ***
## Mujeres      -1.298e-02  2.074e-02  -0.626 0.532507    
## Urbano        4.829e-03  2.546e-03   1.897 0.060352 .  
## ExpectVida    1.283e-04  2.002e-02   0.006 0.994897    
## NeontlMort    2.653e-03  8.600e-03   0.308 0.758299    
## DisMort      -1.046e-02  1.474e-02  -0.710 0.479353    
## Lesion        1.047e-02  1.243e-02   0.842 0.401434    
## EnfTrans     -1.180e-03  8.424e-03  -0.140 0.888823    
## PBI          -4.303e-03  3.634e-03  -1.184 0.238827    
## Tuberculosis -7.273e-04  3.858e-04  -1.885 0.061924 .  
## Diabetes      1.807e-02  1.182e-02   1.529 0.129119    
## Medicos      -3.670e-02  6.256e-02  -0.587 0.558564    
## Camas        -1.134e-02  2.785e-02  -0.407 0.684717    
## ImmunSaramp  -5.581e-03  3.375e-03  -1.653 0.100971    
## TempMarzo    -5.852e-03  5.937e-03  -0.986 0.326348    
## HipTen        3.743e-03  5.562e-03   0.673 0.502362    
## BCG          -1.680e-01  1.718e-01  -0.978 0.330183    
## Tiempo        2.712e-03  2.270e-03   1.195 0.234723    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 0.3624 
## Multiple R-squared:  0.7963, Adjusted R-squared:  0.7555 
## Convergence in 36 IRWLS iterations
plot(modelo_MM_covid_V2$residuals,modelo_MM_covid_V2$rweights, 
     main = "REGRESION ROBUSTA 1", xlab = "RESIDUOS", 
     ylab = "RWEIGHTS",col = "blue")

Vemos que con este método no vemos 2 valores atípicos como en el anterior, sino que vemos 5 de ellos. Desde nuestro punto de vista resultaría muy “costoso” dropear todos estos valores, por lo que no quisieramos recurrir a esta opción. De todas formas vamos a correr un modelo quitando estos datos para analizar si esto mejora la performance del mismo y, de hacerlo, dimensionar en cuanto lo hace.

df6 = df6[(df6$geoId!='US')&(df6$geoId!='CN')&(df6$geoId!='IN')&(df6$geoId!='QA')&(df6$geoId!='JP'),]# DROPEAMOS LOS 5 OUTLIERS QUE DETECTO MM
modelo_l10muertes_V3 = lm(l10muertes.permil~. , data = df6[,-1])
summary(modelo_l10muertes_V3)
## 
## Call:
## lm(formula = l10muertes.permil ~ ., data = df6[, -1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.69790 -0.18409 -0.00631  0.16924  0.83687 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.403e+00  2.606e+00   0.538 0.591442    
## casos.permil  3.127e-04  3.124e-05  10.012  < 2e-16 ***
## Pobla80      -3.494e-02  9.670e-02  -0.361 0.718563    
## Pobla65       6.813e-02  3.126e-02   2.179 0.031443 *  
## PoblaMid     -8.959e-03  1.354e-02  -0.662 0.509491    
## PoblaData     2.931e-01  8.317e-02   3.524 0.000620 ***
## PoblaDens    -1.919e-02  5.326e-03  -3.603 0.000475 ***
## Mujeres      -1.391e-02  1.944e-02  -0.715 0.475819    
## Urbano        3.608e-03  2.314e-03   1.559 0.121912    
## ExpectVida    3.455e-03  1.845e-02   0.187 0.851793    
## NeontlMort    4.801e-03  7.891e-03   0.608 0.544176    
## DisMort      -7.163e-03  1.315e-02  -0.545 0.587176    
## Lesion        1.598e-02  1.145e-02   1.396 0.165653    
## EnfTrans     -2.993e-03  7.875e-03  -0.380 0.704639    
## PBI          -4.247e-03  3.112e-03  -1.365 0.175149    
## Tuberculosis -4.938e-04  3.520e-04  -1.403 0.163498    
## Diabetes      1.302e-02  1.099e-02   1.184 0.238944    
## Medicos      -1.736e-03  5.629e-02  -0.031 0.975455    
## Camas        -2.573e-02  2.531e-02  -1.017 0.311615    
## ImmunSaramp  -4.772e-03  3.106e-03  -1.536 0.127344    
## TempMarzo    -3.039e-03  5.436e-03  -0.559 0.577186    
## HipTen        5.621e-03  5.219e-03   1.077 0.283843    
## BCG          -2.087e-01  1.631e-01  -1.280 0.203301    
## Tiempo        1.904e-03  2.129e-03   0.894 0.373299    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3376 on 110 degrees of freedom
## Multiple R-squared:  0.8364, Adjusted R-squared:  0.8021 
## F-statistic: 24.44 on 23 and 110 DF,  p-value: < 2.2e-16
wise_modelo_muertes_permil_V4 = step(modelo_l10muertes_V3, direction = c("both"),
                                                         trace = F)
formula_elegida = wise_modelo_muertes_permil_V4$call$formula #Nos guardamos las variables seleccionadas para utilizarlas luego
summary(wise_modelo_muertes_permil_V4)
## 
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla65 + PoblaData + 
##     PoblaDens + Urbano + Tuberculosis + Camas + ImmunSaramp + 
##     BCG, data = df6[, -1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.67764 -0.17239  0.00686  0.16236  0.87887 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.795e-01  2.814e-01   3.125 0.002215 ** 
## casos.permil  2.902e-04  2.516e-05  11.532  < 2e-16 ***
## Pobla65       5.165e-02  7.238e-03   7.136 7.10e-11 ***
## PoblaData     2.718e-01  7.021e-02   3.871 0.000174 ***
## PoblaDens    -2.484e-02  4.317e-03  -5.754 6.43e-08 ***
## Urbano        3.467e-03  1.686e-03   2.056 0.041857 *  
## Tuberculosis -7.638e-04  2.510e-04  -3.043 0.002862 ** 
## Camas        -3.255e-02  1.821e-02  -1.787 0.076329 .  
## ImmunSaramp  -5.360e-03  2.492e-03  -2.150 0.033461 *  
## BCG          -2.128e-01  1.486e-01  -1.432 0.154605    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3298 on 124 degrees of freedom
## Multiple R-squared:  0.824,  Adjusted R-squared:  0.8112 
## F-statistic: 64.48 on 9 and 124 DF,  p-value: < 2.2e-16

Vemos que el modelo ajustado por OLS mejora en un 1.6% (0.7978 vs 0.8112). Vamos a volver a ajustar una regresión robusta para inspeccionar nuevamente el gráfico de residuos que disparó este approach, para ello, vamos a usar las mismas variables que surgieron del último proceso de selección de variables.

modelo_MM_covid_V2 = lmrobdetMM(formula_elegida, 
                             data = select_if(df6,is.numeric), control = control)

summary(modelo_MM_covid_V2)
## 
## Call:
## lmrobdetMM(formula = formula_elegida, data = select_if(df6, is.numeric), 
##     control = control)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73780 -0.16801  0.01109  0.18404  0.86496 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.229e-01  3.031e-01   3.045 0.002843 ** 
## casos.permil  3.058e-04  2.967e-05  10.307  < 2e-16 ***
## Pobla65       5.021e-02  8.032e-03   6.251  6.0e-09 ***
## PoblaData     2.806e-01  7.805e-02   3.595 0.000466 ***
## PoblaDens    -2.533e-02  4.468e-03  -5.670  9.5e-08 ***
## Urbano        3.517e-03  1.873e-03   1.877 0.062859 .  
## Tuberculosis -8.000e-04  2.715e-04  -2.947 0.003834 ** 
## Camas        -2.320e-02  2.092e-02  -1.109 0.269644    
## ImmunSaramp  -6.583e-03  2.701e-03  -2.437 0.016209 *  
## BCG          -1.824e-01  1.565e-01  -1.166 0.246015    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 0.3118 
## Multiple R-squared:  0.8243, Adjusted R-squared:  0.8116 
## Convergence in 32 IRWLS iterations
plot(modelo_MM_covid_V2$residuals,modelo_MM_covid_V2$rweights,
     main = "REGRESION ROBUSTA 2", xlab = "RESIDUOS", 
     ylab = "RWEIGHTS",col = "blue")

Vemos que tanto los beta no varían mucho contra el modelo de OLS y en el gráfico de rweights vs residuals tenemos una distribución que no indica presencia de outliers. Este termina siendo nuestro mejor modelo por el momento pero, como desventaja, debemos comprender que dropeamos 5 países (Canadá, Qatar, Estados Unidos, India y Japón) podría tener su costo al momento de intentar extrapolar algún tipo de insight que surja del presente trabajo. De todas formas, aun nos queda un approach por probar, el mismo consiste en buscar alguna transformación de la variable target o las covariables regresoras que nos permita subsanar el problema de outliers e intentar mejorar la performance del modelo.

Prueba 4: transformaciones

Como vimos en puntos anteriores una variable que presenta valores extremos es la de casos.permil, por ello resulta lógico probar sobre ella una transformación del tipo logarítmica pero, para corroborar nuestra hipótesis, vamos a correr un test de Box-Tidwell sobre ella para ver qué tipo de transformación nos sugiere

df7 = COVID[, -c(2, 3, 4, 5, 8, 9, 21, 30, 31, 34, 36)]
car::boxTidwell(l10muertes.permil~casos.permil, data=df7[,-1])
##  MLE of lambda Score Statistic (z)  Pr(>|z|)    
##          0.125             -14.734 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## iterations =  10

Efectivamente, el test Box-Tidwell nos sugiere hacer una transformación logarítmica (dado que un lambda=0 cae dentro de sus intervalos de confianza)

df7$casos.permil = log(df7$casos.permil)
modelo_l10muertes_outliers = lm(l10muertes.permil~. , data = df7[,-1])
summary(modelo_l10muertes_outliers)
## 
## Call:
## lm(formula = l10muertes.permil ~ ., data = df7[, -1])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.78821 -0.19981  0.01797  0.15267  0.71851 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.676e+00  2.280e+00  -1.174 0.242844    
## casos.permil  2.956e-01  2.232e-02  13.244  < 2e-16 ***
## Pobla80       6.935e-02  8.350e-02   0.830 0.407998    
## Pobla65       2.770e-02  2.760e-02   1.003 0.317740    
## PoblaMid      2.612e-03  1.188e-02   0.220 0.826429    
## PoblaData     2.796e-02  1.835e-02   1.524 0.130248    
## PoblaDens    -1.873e-02  4.769e-03  -3.927 0.000147 ***
## Mujeres       3.448e-02  1.597e-02   2.159 0.032948 *  
## Urbano       -5.229e-04  2.087e-03  -0.251 0.802598    
## ExpectVida    5.198e-03  1.599e-02   0.325 0.745660    
## NeontlMort   -4.977e-03  7.136e-03  -0.697 0.486952    
## DisMort       1.696e-03  1.172e-02   0.145 0.885178    
## Lesion        4.963e-03  1.010e-02   0.491 0.624146    
## EnfTrans      6.911e-03  7.071e-03   0.977 0.330411    
## PBI           3.723e-03  2.544e-03   1.463 0.146102    
## Tuberculosis -7.847e-05  3.133e-04  -0.250 0.802682    
## Diabetes     -4.676e-03  9.750e-03  -0.480 0.632407    
## Medicos      -3.004e-02  4.705e-02  -0.639 0.524401    
## Camas        -5.655e-02  2.012e-02  -2.810 0.005825 ** 
## ImmunSaramp  -5.811e-04  2.837e-03  -0.205 0.838076    
## TempMarzo    -2.979e-03  4.790e-03  -0.622 0.535304    
## HipTen       -8.175e-03  4.681e-03  -1.746 0.083439 .  
## BCG          -3.779e-01  1.397e-01  -2.706 0.007848 ** 
## Tiempo        3.096e-03  1.903e-03   1.627 0.106491    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3067 on 115 degrees of freedom
## Multiple R-squared:  0.8632, Adjusted R-squared:  0.8358 
## F-statistic: 31.55 on 23 and 115 DF,  p-value: < 2.2e-16
wise_modelo_muertes_permil_V4 = step(modelo_l10muertes_outliers, direction = c("both"),
                                                         trace = F)
formula_elegida = wise_modelo_muertes_permil_V4$call$formula
summary(wise_modelo_muertes_permil_V4)
## 
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla65 + PoblaData + 
##     PoblaDens + Mujeres + EnfTrans + PBI + Camas + HipTen + BCG + 
##     Tiempo, data = df7[, -1])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7885 -0.1882  0.0072  0.1472  0.7662 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.162721   0.624221  -3.465 0.000724 ***
## casos.permil  0.294724   0.019165  15.378  < 2e-16 ***
## Pobla65       0.048674   0.008718   5.583 1.37e-07 ***
## PoblaData     0.024752   0.016943   1.461 0.146528    
## PoblaDens    -0.019352   0.004019  -4.815 4.10e-06 ***
## Mujeres       0.032079   0.011820   2.714 0.007571 ** 
## EnfTrans      0.004834   0.002032   2.379 0.018839 *  
## PBI           0.004142   0.002167   1.911 0.058218 .  
## Camas        -0.052347   0.015767  -3.320 0.001175 ** 
## HipTen       -0.009433   0.004136  -2.280 0.024251 *  
## BCG          -0.401135   0.125518  -3.196 0.001759 ** 
## Tiempo        0.002724   0.001681   1.621 0.107580    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2956 on 127 degrees of freedom
## Multiple R-squared:  0.8596, Adjusted R-squared:  0.8475 
## F-statistic: 70.69 on 11 and 127 DF,  p-value: < 2.2e-16

Hemos logrado encontrar un mejor modelo, aumentando el r2-ajustado en 4.5% (0.8116 vs 0.8475) y, además, no hemos dropeado ninguna instancia, lo que es súmamente positivo. Vamos a ver los gráfico de residuos para asegurarnos de que no haya valores atipicos.

plot(wise_modelo_muertes_permil_V4)

ks.test(wise_modelo_muertes_permil_V4$residuals, "pnorm", mean(wise_modelo_muertes_permil_V4$residuals),
                                        sd(wise_modelo_muertes_permil_V4$residuals))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  wise_modelo_muertes_permil_V4$residuals
## D = 0.065046, p-value = 0.5988
## alternative hypothesis: two-sided

Pareciera que tenemos distribucion normal y homocedasticidad en los residuos y no tenemos presencia de ningún outlier. Por último, vamos a correr un modelo robusto para comparar los beta con el modelo OLS y para ver el gráfico de residuos vs r-weigths

modelo_MM_covid_log = lmrobdetMM(formula_elegida, 
                             data = select_if(df7,is.numeric), control = control)

summary(modelo_MM_covid_log)
## 
## Call:
## lmrobdetMM(formula = formula_elegida, data = select_if(df7, is.numeric), 
##     control = control)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.79458 -0.17096  0.01575  0.16144  0.85128 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.972730   0.690760  -2.856  0.00502 ** 
## casos.permil  0.275112   0.020656  13.319  < 2e-16 ***
## Pobla65       0.049847   0.009160   5.442 2.61e-07 ***
## PoblaData     0.025107   0.017233   1.457  0.14761    
## PoblaDens    -0.018085   0.004013  -4.506 1.48e-05 ***
## Mujeres       0.028493   0.013021   2.188  0.03048 *  
## EnfTrans      0.005004   0.002146   2.331  0.02132 *  
## PBI           0.004363   0.002297   1.900  0.05975 .  
## Camas        -0.050766   0.016571  -3.064  0.00267 ** 
## HipTen       -0.008630   0.004282  -2.016  0.04594 *  
## BCG          -0.448764   0.133414  -3.364  0.00102 ** 
## Tiempo        0.003983   0.001782   2.235  0.02714 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 0.2881 
## Multiple R-squared:  0.8604, Adjusted R-squared:  0.8483 
## Convergence in 35 IRWLS iterations
plot(modelo_MM_covid_log$residuals,modelo_MM_covid_log$rweights, 
     main = "REGRESION ROBUSTA 3", xlab = "RESIDUOS", 
     ylab = "RWEIGHTS",col = "blue")

Claramente no se encuentra presencia de outliers y los beta estimados no presentan cambios significativos, por lo que podemos concluír en que nuestro modelo no se encuentra afectado por outliers significativos. Un punto más en cuanto a la selección del modelo, vemos que el método robusto mejora el r2-ajustado y el RSS vs la versión de OLS pero, como sabemos, este tipo de método tiene intervalos de confianza mayores, por lo que (al ser marginal la mejora), perferimos quedarnos con el modelo OLS.

final_model = wise_modelo_muertes_permil_V4

summary(lm.covid)
## 
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + 
##     Tuberculosis + Camas + TempMarzo + PBI, data = COVID)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.08088 -0.29012 -0.03314  0.27929  1.29267 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.5981046  0.2012630   2.972  0.00352 ** 
## PoblaDens    -0.0111873  0.0063386  -1.765  0.07990 .  
## Pobla80       0.2040913  0.0338369   6.032 1.55e-08 ***
## Urbano        0.0082702  0.0025767   3.210  0.00167 ** 
## Tuberculosis -0.0005815  0.0003431  -1.695  0.09246 .  
## Camas        -0.1101291  0.0263087  -4.186 5.17e-05 ***
## TempMarzo    -0.0143117  0.0057068  -2.508  0.01337 *  
## PBI           0.0062466  0.0027181   2.298  0.02314 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4765 on 131 degrees of freedom
## Multiple R-squared:  0.6239, Adjusted R-squared:  0.6038 
## F-statistic: 31.04 on 7 and 131 DF,  p-value: < 2.2e-16
summary(final_model)
## 
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla65 + PoblaData + 
##     PoblaDens + Mujeres + EnfTrans + PBI + Camas + HipTen + BCG + 
##     Tiempo, data = df7[, -1])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7885 -0.1882  0.0072  0.1472  0.7662 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.162721   0.624221  -3.465 0.000724 ***
## casos.permil  0.294724   0.019165  15.378  < 2e-16 ***
## Pobla65       0.048674   0.008718   5.583 1.37e-07 ***
## PoblaData     0.024752   0.016943   1.461 0.146528    
## PoblaDens    -0.019352   0.004019  -4.815 4.10e-06 ***
## Mujeres       0.032079   0.011820   2.714 0.007571 ** 
## EnfTrans      0.004834   0.002032   2.379 0.018839 *  
## PBI           0.004142   0.002167   1.911 0.058218 .  
## Camas        -0.052347   0.015767  -3.320 0.001175 ** 
## HipTen       -0.009433   0.004136  -2.280 0.024251 *  
## BCG          -0.401135   0.125518  -3.196 0.001759 ** 
## Tiempo        0.002724   0.001681   1.621 0.107580    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2956 on 127 degrees of freedom
## Multiple R-squared:  0.8596, Adjusted R-squared:  0.8475 
## F-statistic: 70.69 on 11 and 127 DF,  p-value: < 2.2e-16
plot(lm.covid)

plot(final_model)

Performance: vemos el salto significativo que hay entre los 2 modelos, logrando en el segundo un incremento en la performance marcado por: > un aumento del r2-ajustado de 0.6038 a 0.8475 (40%) > un descenso del RSS de 0.4765 a 0.2956 (62%) > la eliminación de la multicolinealidad > el manejo de outliers.

Cambios de variables: es destacable el cambio de variables seleccionadas y betas estimados, en el segundo modelo vemos que casos.permil es una variable preponderante (por ser la variable más grande y tener a su vez el beta estimado más grande) y en el primer modelo no la habíamos incluído. Algo parecido ocurre con la variable BCG, que toma reelevancia en el segundo modelo y en el primero no estaba.

par(mfrow=c(1, 2))
int_prediccion = predict(lm.covid, newdata = COVID, interval = c("predict"))
int_prediccion_2 = predict(wise_modelo_muertes_permil_V4, newdata = df7, interval = c("predict"))

plot(y = lm.covid$fitted.values, 
     x = COVID$l10muertes.permil, 
     col = "darkblue", ylim = c(-2, 5), cex = .8,
     xlab = "VALORES OBSERVADOS", ylab = "VALORES PREDICHOS", main='Intervalos de predicción modelo 1')
# Intervalos de prediccion
points(x = COVID$l10muertes.permil, 
      y = int_prediccion[, 2], col = "darkred", lwd =3, cex = 0.5)
points(x = COVID$l10muertes.permil, 
       y = int_prediccion[, 3], col = "darkred", lwd =3, cex = 0.5)
legend("topleft", c("INTERVALO DE PREDICCION", "OBSERVACION"), fill = c("darkred", "darkblue"))


plot(y = wise_modelo_muertes_permil_V4$fitted.values, 
     x = df7$l10muertes.permil, 
     col = "darkblue", ylim = c(-2, 5), cex = .8,
     xlab = "VALORES OBSERVADOS", ylab = "VALORES PREDICHOS", main='Intervalos de predicción modelo 2')
# Intervalos de prediccion
points(x = df7$l10muertes.permil, 
      y = int_prediccion_2[, 2], col = "darkred", lwd =3, cex = 0.5)
points(x = df7$l10muertes.permil, 
       y = int_prediccion_2[, 3], col = "darkred", lwd =3, cex = 0.5)
legend("topleft", c("INTERVALO DE PREDICCION", "OBSERVACION"), fill = c("darkred", "darkblue"))

En cuanto a los intervalos de predicción, vemos como en nuestro último modelo estos de vuelven más estrechos, acercándose a los valores reales de nuestra variable respuesta.

plot(final_model)


Observamos que un factor muy ponderante para disminuir las muertes debidas al COVID19, fue la estrategia de INMUNIZACION, por lo que sugeririamos a las autoridades que opten por una vacunacion TOTAL y no sectorizada, y se difundan mas publicidades para concientizar sobre los beneficios de vacunarse. Tambien aumentar la cantidad de campañas y centros de vacunacion; incluso de ser posible, dar incentivos a los que aun conservan dudas sobre si vacunarse, para que lo hagan.

Otro factor que influyo en gran medida fue la cantidad de casos por cada mil habitantes, por lo que recomendamos avanzar con estrategias de aislamiento preventivo y campañas de concientización de los cuidados escenciales como el uso de barbijo, higiene correcta de manos y distancia social.

Por último, tambien observamos que la poblacion mas afectada por el virus, fue la de mayor edad (mas de 65 años). Sugerimos implementar mas politicas de concientizacion a la sociedad sobre la vulneravilidad de este sector, para que se extremen más los recaudos. Tambien sugerimos agrandar los programas publicos de asistencia para personas mayores y lanzar voluntariados a la sociedad en general, para suplir esta demanda. De esta forma, la persona mayor puede disminuir la chance de contagio, evitando ir a lugares con mucha densidad de gente (supermercados, bancos, farmacias, etc), para cubrir sus necesidades basicas.